R Notebook

This is an R Markdown Notebook. When you execute code within the notebook, the results appear beneath the code.

Try executing this chunk by clicking the Run button within the chunk or by placing your cursor inside it and pressing Ctrl+Shift+Enter.

getwd()

[1] "C:/Users/rchaz/Downloads"

Add a new chunk by clicking the Insert Chunk button on the toolbar or by pressing Ctrl+Alt+I.

When you save the notebook, an HTML file containing the code and output will be saved alongside it (click the Preview button or press Ctrl+Shift+K to preview the HTML file).

The preview shows you a rendered HTML copy of the contents of the editor. Consequently, unlike Knit, Preview does not run any R code chunks. Instead, the output of the chunk when it was last run in the editor is displayed.

#analyse du jeu de données MP graminées
#en R cette fois

install.packages("dplyr")

WARNING: Rtools is required to build R packages but is not currently installed. Please download and install the appropriate version of Rtools before proceeding:

https://cran.rstudio.com/bin/windows/Rtools/
Installing package into ‘C:/Users/rchaz/AppData/Local/R/win-library/4.3’
(as ‘lib’ is unspecified)
trying URL 'https://cran.rstudio.com/bin/windows/contrib/4.3/dplyr_1.1.2.zip'
Content type 'application/zip' length 1551420 bytes (1.5 MB)
downloaded 1.5 MB

package ‘dplyr’ successfully unpacked and MD5 sums checked

The downloaded binary packages are in
    C:\Users\rchaz\AppData\Local\Temp\RtmpAHim9b\downloaded_packages

library(rlang)
library(dplyr)

Attaching package: ‘dplyr’

The following object is masked from ‘package:signal’:

    filter

The following objects are masked from ‘package:stats’:

    filter, lag

The following objects are masked from ‘package:base’:

    intersect, setdiff, setequal, union

#installation des librairies
install.packages("readxl")

WARNING: Rtools is required to build R packages but is not currently installed. Please download and install the appropriate version of Rtools before proceeding:

https://cran.rstudio.com/bin/windows/Rtools/
Installing package into ‘C:/Users/rchaz/AppData/Local/R/win-library/4.3’
(as ‘lib’ is unspecified)
trying URL 'https://cran.rstudio.com/bin/windows/contrib/4.3/readxl_1.4.3.zip'
Content type 'application/zip' length 1194912 bytes (1.1 MB)
downloaded 1.1 MB

package ‘readxl’ successfully unpacked and MD5 sums checked

The downloaded binary packages are in
    C:\Users\rchaz\AppData\Local\Temp\RtmpAHim9b\downloaded_packages

install.packages("ggplot2")

WARNING: Rtools is required to build R packages but is not currently installed. Please download and install the appropriate version of Rtools before proceeding:

https://cran.rstudio.com/bin/windows/Rtools/
Installing package into ‘C:/Users/rchaz/AppData/Local/R/win-library/4.3’
(as ‘lib’ is unspecified)
trying URL 'https://cran.rstudio.com/bin/windows/contrib/4.3/ggplot2_3.4.3.zip'
Content type 'application/zip' length 3329305 bytes (3.2 MB)
downloaded 3.2 MB

package ‘ggplot2’ successfully unpacked and MD5 sums checked

The downloaded binary packages are in
    C:\Users\rchaz\AppData\Local\Temp\RtmpAHim9b\downloaded_packages

install.packages("caret")

WARNING: Rtools is required to build R packages but is not currently installed. Please download and install the appropriate version of Rtools before proceeding:

https://cran.rstudio.com/bin/windows/Rtools/
Installing package into ‘C:/Users/rchaz/AppData/Local/R/win-library/4.3’
(as ‘lib’ is unspecified)
trying URL 'https://cran.rstudio.com/bin/windows/contrib/4.3/caret_6.0-94.zip'
Content type 'application/zip' length 3576720 bytes (3.4 MB)
downloaded 3.4 MB

package ‘caret’ successfully unpacked and MD5 sums checked

The downloaded binary packages are in
    C:\Users\rchaz\AppData\Local\Temp\RtmpAHim9b\downloaded_packages

install.packages("MASS") 

Error in install.packages : Updating loaded packages

install.packages("reshape2") 

WARNING: Rtools is required to build R packages but is not currently installed. Please download and install the appropriate version of Rtools before proceeding:

https://cran.rstudio.com/bin/windows/Rtools/
Installing package into ‘C:/Users/rchaz/AppData/Local/R/win-library/4.3’
(as ‘lib’ is unspecified)
trying URL 'https://cran.rstudio.com/bin/windows/contrib/4.3/reshape2_1.4.4.zip'
Content type 'application/zip' length 454253 bytes (443 KB)
downloaded 443 KB

package ‘reshape2’ successfully unpacked and MD5 sums checked

The downloaded binary packages are in
    C:\Users\rchaz\AppData\Local\Temp\RtmpAHim9b\downloaded_packages

install.packages("reshape") 

WARNING: Rtools is required to build R packages but is not currently installed. Please download and install the appropriate version of Rtools before proceeding:

https://cran.rstudio.com/bin/windows/Rtools/
Installing package into ‘C:/Users/rchaz/AppData/Local/R/win-library/4.3’
(as ‘lib’ is unspecified)
trying URL 'https://cran.rstudio.com/bin/windows/contrib/4.3/reshape_0.8.9.zip'
Content type 'application/zip' length 170892 bytes (166 KB)
downloaded 166 KB

package ‘reshape’ successfully unpacked and MD5 sums checked

The downloaded binary packages are in
    C:\Users\rchaz\AppData\Local\Temp\RtmpAHim9b\downloaded_packages

install.packages("signal") 

Error in install.packages : Updating loaded packages

#Load required libraries
library(ggplot2)
library(readxl)
library(caret)

Loading required package: lattice

Attaching package: ‘caret’

The following object is masked from ‘package:pls’:

    R2

library(signal)
library(pls)

# Create an empty data frame
df <- data.frame()

# Read the Excel data
data <- read_excel('MPall2000600cm.xlsx')

# Display the first 50 rows of data
head(data, 50)

#49 échantillons... c'st tout ce que j'ai pu retrouver.
#et encore : il y a 10 lignines dedans. Vu qu'on se concentre sur les MP, ça réduit le br d'échantillons à 38...
#on ne peut pas faire une prédiction dans les règles avec si peu. Je ne sais pas si je pourrai même obtenir des performances satisfaisantes ave si peu.

#Ouliers attendus : 9F, FF1, FF3  (6, 7, 8)

# Analyze variable types
variable_types <- table(sapply(data, class))

# Plot pie chart for variable types
pie(variable_types, labels = names(variable_types), main = "Variable Types")

# Get the counts of variable types
variable_counts <- table(sapply(data, class))

# Print the counts of variable types
print(variable_counts)

character   numeric 
        1       737 

library(reshape2)

# Drop rows with indices 6, 7, and 8
data <- data[-c(6, 7, 8), ]

# Display the first 40 rows of data
head(data, 40)

# Drop row index 6
data <- data[-6, ]

# Select columns for X
X <- data[, 12:ncol(data)]

# Select columns for Y in PLS1
Y_pcoum <- data[, 10]
Y_SG <- data[, 7]
Y_syr <- data[, 8]
Y_gua <- data[, 9]
Y_oses <- data[, 2]
Y_LK <- data[, 6]

# Select columns for Y in PLS2
Y <- data[, c(2, 6, 7, 10)]

# Display the first 50 rows of X dataset
head(X, 50)

# Display the first few rows of Y_pcoum column
head(Y, 50)

# Create a heatmap to visualize missing values
install.packages("ggplot2")

Error in install.packages : Updating loaded packages

library(ggplot2)

missing_values <- is.na(data)
missing_values_melted <- melt(missing_values)

ggplot(data = missing_values_melted, aes(x = Var2, y = Var1, fill = value)) +
  geom_tile() +
  scale_fill_manual(values = c("TRUE" = "red", "FALSE" = "white")) +labs(title = "Missing Values Heatmap")

# Drop rows with indices 6, 7, and 8
#data <- data[-c(6, 7, 8), ]

# Display the first 40 rows of data
head(data, 40)

# Drop row index 6
#data <- data[-6, ]

# Select columns for X
X <- data[, 12:ncol(data)]

# Select columns for Y in PLS1
Y_pcoum <- data[, 10]
Y_SG <- data[, 7]
Y_syr <- data[, 8]
Y_gua <- data[, 9]
Y_oses <- data[, 2]
Y_LK <- data[, 6]

# Select columns for Y in PLS2
Y <- data[, c(2, 6, 7, 10)]

# Display the first 50 rows of X dataset
head(X, 50)

# Display the first few rows of Y_pcoum column
head(Y_pcoum, 50)

# Display the first few rows of Y dataset
head(Y)

#proChaine étape : le prétraitement des spectres


# Create a vector for wavelength values
wl <- seq(0, 727, 1)

# Plot the transposed data matrix
matplot(X, type = "l", xlab = "Wavelength", ylab = "Spectral Value", 
        col = 1:ncol(X), lty = 1)

dim(Y)

[1] 35  4

dim(X)

[1]  35 727

# Define the SNV normalization function
snv <- function(input_data) {
  output_data <- matrix(0, nrow = nrow(input_data), ncol = ncol(input_data))
  
  for (i in 1:nrow(input_data)) {
    output_data[i,] <- (input_data[i,] - mean(input_data[i,])) / sd(input_data[i,])
  }
  
  return(output_data)
}

# Convert X to a matrix
X <- as.matrix(X)

# Apply SNV normalization
Xsnv <- snv(X)

# Plot the normalized spectral values
matplot(t(Xsnv), type = "l", xlab = "Wavelength", ylab = "Normalized Spectral Value", 
        col = 1:ncol(Xsnv), lty = 1)

dim(Xsnv)

[1]  35 727

# Load the required library
install.packages("signal") 

Error in install.packages : Updating loaded packages

library(signal)

# Define Savitzky-Golay parameters
window_size <- 17
poly_order <- 2

# Transpose the matrix to work on columns
X_transposed <- t(X)

# Apply the Savitzky-Golay filter for second derivative to columns
X2_transposed <- apply(X_transposed, 2, function(col) {
  deriv_filter <- sgolayfilt(c(rep(0, 8), -2, rep(0, 8)), p = 2, n = 17)
  filtered_col <- convolve(col, deriv_filter, type = "filter")
  return(filtered_col)
})

# Transpose the result back to the original orientation
X2 <- t(X2_transposed)

# Plot the second derivative values
matplot(t(X2), type = "l", xlab = "Wavelength", ylab = "Second Derivative Value", 
        col = 1:ncol(X2), lty = 1)

dim(X2)

[1]  35 711

# Apply the Savitzky-Golay filter for second derivative
X2snv <- snv(X2)

dim(X2snv)

[1]  35 711

#tentons une PLS
library(pls)


# Initialize empty vectors to store results
R2_cal <- numeric()
R2_cv <- numeric()
RMSE_cal <- numeric()
RMSE_cv <- numeric()
Var_cal <- numeric()

Components <- 1:20

# Combine the predictor variables and the response variable into a single data frame
data_df <- data.frame(Y_pcoum = unlist(Y_pcoum), Xsnv)


# Calibration
PLS_pcoum <- plsr(Y_pcoum ~ ., ncomp = 7, data = data_df)
Ypred <- predict(PLS_pcoum, newdata = data_df)
  
# Extract the appropriate dimension of Ypred
Ypred_flat <- Ypred[, 1, ]
  
#reprenons nos calculs
R2 <- cor(Ypred_flat, data_df$Y_pcoum)^2
RMSE <- sqrt(mean((Ypred_flat - data_df$Y_pcoum)^2))
 
#résultats 
cat("Results for 7 components\n")

Results for 7 components

cat("R2:", R2, "\n")

R2: 0.5911774 0.6985943 0.854471 0.9495035 0.9688896 0.9771441 0.9819514 

cat("RMSE:", RMSE, "\n\n")

RMSE: 3.036214 

# Cross-validation
Y_cv <- predict(PLS_pcoum, newdata = data.frame(Xsnv), validation = "CV")

# Extract the appropriate dimension of Ypred
Y_cv_flat <- Y_cv[, 1, ]
  
# Calculate scores
R2_cv2 <- cor(Y_cv_flat, data_df$Y_pcoum)^2
RMSE_cv <- sqrt(mean((Y_cv_flat - data_df$Y_pcoum)^2))
  
# Print the results
cat("Results for 7 components\n")

Results for 7 components

cat("R2_cv2:", R2_cv2, "\n")

R2_cv2: 0.5911774 0.6985943 0.854471 0.9495035 0.9688896 0.9771441 0.9819514 

cat("RMSE_cv:", RMSE_cv, "\n\n")

RMSE_cv: 3.036214 

#tentons une autre aporoche pour la cross validation.

# Calibration
PLS_pcoum <- plsr(Y_pcoum ~ ., ncomp = 7, data = data_df)
Ypred <- predict(PLS_pcoum, newdata = data_df)

# Extract the appropriate dimension of Ypred
Ypred_flat <- Ypred[, 1, ]

# Perform manual cross-validation
set.seed(123)  # For reproducibility
num_samples <- nrow(data_df)
num_folds <- 10
fold_size <- num_samples / num_folds

# Initialize vectors to store cross-validation results
R2_cv <- numeric()
RMSE_cv <- numeric()

for (fold in 1:num_folds) {
  # Divide data into training and validation sets
  val_start <- ((fold - 1) * fold_size) + 1
  val_end <- fold * fold_size
  val_indices <- val_start:val_end
  train_indices <- setdiff(1:num_samples, val_indices)
  
  # Fit PLS model on training data
  PLS_pcoum_fold <- plsr(Y_pcoum ~ ., ncomp = 7, data = data_df[train_indices, ])
  
  # Predict on validation data
  Y_cv_fold <- predict(PLS_pcoum_fold, newdata = data_df[val_indices, ])
  
  # Extract the appropriate dimension of Y_cv_fold
  Y_cv_flat <- Y_cv_fold[, 1, ]
  
  # Calculate scores for this fold
  R2_cv_fold <- cor(Y_cv_flat, data_df$Y_pcoum[val_indices])^2
  RMSE_cv_fold <- sqrt(mean((Y_cv_flat - data_df$Y_pcoum[val_indices])^2))
  
  # Store the results for this fold
  R2_cv <- c(R2_cv, R2_cv_fold)
  RMSE_cv <- c(RMSE_cv, RMSE_cv_fold)
}


# Calculate average R2 for cross-validation
average_R2_cv <- mean(R2_cv)

# Calculate average RMSE for cross-validation
average_RMSE_cv <- mean(RMSE_cv)

# Print the cross-validation results
cat("Cross-validation results:\n")

Cross-validation results:

cat("R2_cv:", average_R2_cv, "\n")

R2_cv: 0.7519379 

cat("RMSE_cv:", average_RMSE_cv, "\n")

RMSE_cv: 2.781517 

#j'arrete le script ici. Mon objectif est de monter que je peux traduire un script python en R.
#Le reste du script étant de la PLS sur d'autres paramètres et avec les spectres dérivés, c'est assez répétitif.
#